function [ A ] = assemblyMatrix( mesh )
%ASSEMBLYMATRIX Summary of this function goes here
%   Detailed explanation goes here

nTri = size(mesh.f, 1);

A = diag(ones(3*nTri, 1)) * 0.5;

for i = 1 : nTri
    
   p0 = mesh.v(mesh.f(i, 1), :);
   p1 = mesh.v(mesh.f(i, 2), :);
   p2 = mesh.v(mesh.f(i, 3), :);
   
   p = (p0 + p1 + p2) / 3.0;
   
   for j = 1 : nTri
       
       q0 = mesh.v(mesh.f(j, 1), :);
       q1 = mesh.v(mesh.f(j, 2), :);
       q2 = mesh.v(mesh.f(j, 3), :);
       
       T = triarea(q0, q1, q2);
              
       qn = mesh.fn(j, :);
       
       v0 = evaluateGradientQuadraturePoint(p, q0, qn);
       v1 = evaluateGradientQuadraturePoint(p, q1, qn);
       v2 = evaluateGradientQuadraturePoint(p, q2, qn);
                     
       A(3*i-2:3*i, 3*j-2:3*j) = A(3*i-2:3*i, 3*j-2:3*j) + (v0 + v1 + v2) * T / 3.0;
       
   end
end

end

